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Abstract 

We present here a short and subjective review of methods for calculating charge trans¬ 
fer couplings. Although we mostly focus on Density Functional Theory, we discuss a small 
subset of semiempirical methods as well as the adiabatic-to-diabatic transformation methods 
typically coupled with wavefunction-based electronic structure calculations. In this work, we 
will present the reader with a critical assessment of the regimes that can be modelled by the 
various methods - their strengths and weaknesses. In order to give a feeling about the practical 
aspects of the calculations, we also provide the reader with a practical protocol for running 
coupling calculations with the recently developed FDE-ET method. 
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1 Introduction 


Charge transfer (CT) between molecular species play vital roles in processes that occur in biology 
such as protein communication, 1-5 respiratory systems in the mitochondria, 6 oxidative damage on 
DNA, 7 9 photo synthetic cycles, 10,11 as well as in materials science conduction in organic semi¬ 
conductors. 12,13 In order to achieve an accurate modeling of these processes in the simulations, 
one needs to include several levels of complexity, which in most instances lead to considering 
model systems featuring hundreds of atoms and an even larger number of electrons. The large 
system sizes preclude the use of high-level wavefunction-based quantum-chemical methods. For 
this reason, researchers worldwide have invested a great deal of effort in developing approximate, 
fast, yet still accurate methods for describing CT reactions. Methods based on Density-Functional 
Theory have in recent years become competitive in regards to the accuracy while still maintaining 
a generally low computational cost. 

Marcus theory 14,15 is perhaps the most applicable theory for modeling a CT process. This 
theory was originally derived under three main approximations. First, a CT event is thought of 
in terms of a two-dimensional basis set (donor and acceptor). The interaction matrix element of 
the Hamiltonian is the central quantity in determining the probability of a transition in populations 
from the basis function representing the donor state to one representing the acceptor state. This in¬ 
teraction is the electronic coupling Vqa of the two electronic states involved in the CT reaction. 16,17 
Second, it relies on the Condon approximation, 18,19 in which the electronic coupling is considered 
to be independent of the nuclear motion when the transfer occurs. Third, reactants and products 
are modeled as being enclosed by spheres on which the polarization of the solvent is represented 
as a dielectric continuum. 15,2 °" 22 Marcus theory can be summarized as: 23 


(AG+A) 2 

2k n e uk b t 

k ^ = T^ 2 7W- 


(1) 


where A is the reorganization energy, and Vda is the electronic coupling. 

States that most resemble the initial and final states of electron transfer have been often referred 
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to as “diabatic states” 24,25 and their corresponding wavefunctions “diabats”. Although it is known 
that diabatic states have a formal definition, 26,27 it was shown 28 that charge-localized states satisfy 
the requirements for diabatic states for condensed phase electron transfer reactions. 

Several approaches are available in the literature to generate and evaluate Hamiltonian matrix 
elements with wavefunctions of charge-localized, diabatic states. They differ in the level of the¬ 
ory used in the calculation and in the way localized electronic structures are created. 15 > 2 5,26,29-3i 
When wavefunction-based quantum-chemical methods are employed, the framework of the gen¬ 
eralized Mulliken-Hush method (GMH), 29,32-34 is particularly successful. So far, it has been 
used in conjunction with accurate electronic structure methods for small and medium sized sys¬ 
tems. 35-37 As an alternative to GMH and other derived methods, 38,39 additional methods have been 
explored for their applicability in larger systems such as constrained density functional method 
(CDFT), 25,37,40,41 and fragmentation approaches, 42-47 which also include the frozen density em¬ 
bedding (FDE) method. 48,49 

So far, we have mentioned methods that produce all-electron diabatic wavefunctions and cor¬ 
responding Hamiltonian matrix elements. There are two other classes of methods which simplify 
the quantum problem by focusing on the wavefunction of the transferred charge: such as methods 
making use of the frozen core approximation Fragment Orbital methods (FO), and methods that 
assume the charge to be localized on single atomic orbitals. 50 In this work, we will also treat these 
computationally low-cost methods. 

As our group is involved in the development of the Frozen Density Embedding (FDE) for¬ 
mulation of subsystem DFT, this chapter will pay particular attention to the FDE methodology. 
We believe FDE to be a very promising method capable of achieving a good description of the 
electronic coupling in CT reactions, while maintaining a low computational complexity. 

This chapter is divided in two parts: the first part is devoted to the FDE method as well as other 
DFT-based alternatives. The second part covers more accurate methods (wavefunction-based). 
In each of the two parts, we discuss the numerical stability and accuracy of the methods in the 
generation of diabatic states with the overarching goal of obtaining reliable electronic couplings 
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with a contained computational effort. 

We will start with a description of FDE and its ability to generate diabats and to compute 
Hamiltonian matrix elements - the FDE-ET method (ET stands for Electron Transfer). In the 
subsequent section, we will present specific examples of FDE-ET computations to provide the 
reader with a comprehensive view of the performance and applicability of FDE-ET. After FDE has 
been treated, four additional methods to generate diabatic states are presented in order of accuracy: 
CDFT, FODFT, AOM, and Pathways. In order to output a comprehensive presentation, we also 
describe those methods in which wavefunctions methods can be used, in particular GMH and 
other adiabatic-to-diabatic diabatization methods. Finally, we provide the reader with a “protocol” 
for running FDE-ET calculations with the only available implementation of the method in the 
Amsterdam Density Functional software. 51 In closing, we outline our concluding remarks and our 
vision of what the future holds for the field of computational chemistry applyed to electron transfer. 

2 DFT Based Methods 

2.1 The Frozen Density Embedding formalism 

The frozen density-embedding (FDE) formalism 52 developed by Wesolowski and Warshel 52-54 
has been applied to a plethora of chemical problems, for instance, solvent effects on different types 
of spectroscopy, 55-57 magnetic properties, 58-62 excited states, 55,63-66 charge transfer states. 49,67,68 
Computationally, FDE is available for molecular systems in ADF, 51,69 Dalton, 70,71 Q-Chem, 72,73 
and Turbomole 74-76 packages, as well as for molecular periodic systems in CP2K 77,78 and fully 
periodic systems (although in different flavors) in CASTEP, 79,80 Quantum Espresso, 81-83 and 
Abinit. 84,85 

FDE prescribes that the total electron density should be expressed as the sum of subsystem 
electron densities, 53,86-89 this is based on the idea that a molecular system can be more easily 
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approached if it is subdivided into many smaller subsystems. Namely, 


#ofsubsystems 

Plot ( r ) = £ P/(r) 


( 2 ) 


/=! 


As in regular DFT calculations, the electron density of each subsystem is computed by solving 
selfconsistently a Kohn-Sham (KS) like equation per subsystem. These KS like equations read as: 


-V 2 


+ v£s( r ) + vLb( r ) 


0(0/( r ) = e (/)/( r )^)/( 1 


(3) 


Where 0(,-)/(r), £(,•)/ are the molecular orbitals and orbital energies of subsystem /. In Eq.(Eq. (3)) 
we have augmented the Kohn-Sham single particle Hamiltonian by an embedding potential, V ! emh , 
in which are encoded the interactions with the other subsystems. In the following, v 1 , (r) is the 
embedding potential acting on subsystem /: 
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/ 
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(4) 


In the above, T s , E xc and Z a are kinetic and exchange-correlation energy functionals, and the 
nuclear charge, respectively, and Ns is the total number of subsystems considered. In practical FDE 
calculations, the kinetic energy is calculated in terms of orbital free semilocal functionals. This 
approximation is ultimately the biggest difference between an FDE and a full KS-DFT calculation 
of the supersystem. 90-93 As a consequence, the embedding potential becomes inaccurate when 
the subsystems feature a large overlap between their electron densities 83,94,95 (this is because the 
larger the density overlap is, the larger the magnitude of the nonadditive potentials become). In 
FDE, the subsystem KS equations are left to converge to selfconsistency with respect to each other. 
This is often achieved by employing the so-called freeze-and-thaw procedure 69,96 (as done in ADF 
and other molecular codes) or via updating the embedding potential at every SCF cycle as done 
in CP2K 78,97 and Quantum-Espresso. 81,83,98 It is worth noting that FDE scales linearly with the 
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number of subsystems provided that linear scaling methods for the solution of the electrostatic 
problem are employed. 69 

The earliest example of diabatization by FDE was given in Ref. 48 This is shown in Figure 
Figure 1, where the spin densities for a pair of guanines are calculated. KS-DFT of the supersystem 
carried out with semilocal XC functionals fails in the prediction of the spin density. This is because 
the self-interaction error makes the spin density spread on both guanines against the prediction 
given by more accurate theoretical work 35 and experimental studies. 99-101 On the contrary, FDE 
locali z es the charge on a guanine of choice. 

The fact that FDE was able to provide subsystem-localized electronic structures was known 
since its early application to systems with unpaired electrons. 60-62,102 Eater, this ability of FDE 
was explored for the computation of diabatic states for electron transfer 48,49,67,103 and to compute 
hyperfine coupling constants. 60,61 



(a) (b) (c) 


Figure 1: Spin densities of a guanine-cytosine dimer radical cation, (GC)^• (a) KS-DFT 

supramolecular calculation using PW91 functional, (b) FDE calculation considering two subsys¬ 
tems where the left side subsystems (blue contour) is positively charged and (c) FDE calculation 
for four subsystems with one subsystems (blue contour) is positively charged. The nucleobases 
structures and spin densities were taken from ref 48. 


The question that one can raise is why FDE calculations yield charge localized states? We 


7 




provide here four reasons. 


49,103 

1. Orthogonality is not imposed between the molecular orbitals belonging to different subsys¬ 
tems. 

2. FDE calculations can be executed in the monomer basis set. This is known as FDE(m) 
method. 104 

3. FDE calculations are always initiated with a subsystem localized initial guess of the electron 
density. 

4. Electrons of a subsystem, remain localized also because there are repulsive walls in the 
region of the surrounding (frozen) fragments. 

The first reason, is important because it directly removes a bias towards delocalization which 
results due to orthonormalization of the molecular orbitals, as already noted by Dulak and Wesolowski. 
The second and third reasons come together, the lack of basis functions on the surrounding subsys¬ 
tems, does not allow substantial charge transfer between the subsystems. As a consequence, the 
SCF is biased to converge to localized electronic structures. 

The fourth reason makes reference to the approximate nature of the term (also 

known as nonadditive kinetic energy potential which is part of the embedding potential) in the 
region of the frozen fragments (e.g. in the region where pj with 7 ^ / is larger than any other 
subsystem electron density). Approximate nonadditive kinetic energy potentials fail in canceling 
out the attractive potential due to the nuclear charge in the vicinity of the nucleus of the sur¬ 
rounding frozen subsystems, 104,106 and they do not reproduce the exact potential at intermediate 
regions, 1 07_109 especially in the vicinity of an atomic shell. 1 ° 6,109 In that region they cross the ex¬ 
act potential and large potential walls arise. A simplified depiction of this effect is devised in figure 
Figure 2. 

In this scenario, diabatic states can be generated with FDE by performing at least two simula¬ 
tions, one featuring a hole/electron on the donor while the acceptor is neutral and one calculation in 
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Figure 2: Exemplification of the embedding potential at the atomic shells of the surrounding 
subsystems. Figure from reference 103. 


which the charge hole/electron is on the acceptor. The result is two charge localized states, whose, 
densities and Kohn-Sham orbitals are used in a later step in order to build the diabatic Hamiltonian 
and overlap matrices, needed to compute the diabatic coupling matrix element. 

2.1.1 FDE-ET method 

FDE-ET is a methodology which computes Hamiltonian couplings from diabatic states generated 
by an FDE calculation. Electron transfer reaction are usually described in the basis of a two-state 
formalism, 30 taking as basis set two broken-symmetry charge-localized states. This methodology 
can also approach models for the superexchange mechanism, 7,23,110-1 13 where the transfer is still 
modelled by a Two-dimensional basis set but the coupling includes the effect of non-resonant 
bridges states. Figure Figure 3, illustrates the difference between tunneling through the vacuum 
and through a set of bridge states. The bridge could be comprised of one or more molecules, a 
covalent bond or any other type o potential barrier as long as its height is lower than the one when 
vacuum separates donor and acceptor. As it is shown in Figure Figure 3, the higher the potential 
barrier the faster the coupling decays with respect to the donor-acceptor distance. 

In FDE-ET we seek a method capable of computing the Hamiltonian matrix in the basis of 
charge-localized states generated with FDE. First, we have to define the needed matrix elements. 
As diabatic states are not the eigenfunctions of the molecular Hamiltonian, the off-diagonal ele- 
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Energy (eV) 





Figure 3: Energy dependence of a charge transfer process. The off-diagonal element (Dirac 
notation) will decay as the potential well that the charge has to overcome increases. Two cases: for 
vacuum as a potential well we have a faster decay, and when molecules act as a bridge the transfer 
will decay slower. 
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ments of such Hamiltonian are not zero and can be approximated by the following formula 
if Yd and Ya are slater determinants representing the donor and acceptor diabats: 


114,115 


Hda — (Vd\H\Ya) — SdaE 


.(DA) 


(r) 


(5) 


Here H is the molecular electronic Hamiltonian, and p il)A) (r') is the transition density defined as 
p( DA )(r) — (YdI'LIL] <5 (it — *)\Va), with n e being the total number of electrons in the system (i.e. 
the sum of the electron number of all subsystems) and E p {DA ^ (r) is an energy density functional. 
The donor-acceptor overlap matrix elements are found by computing the following determinant: 


Sda = det 


S (DA) 


( 6 ) 


where S ! k >A = ( /J ' | A *) is the transition overlap matrix in terms of the occupied orbitals 

Thus, the transition density is now written in the basis of all occupied orbitals which make up the 
diabatic states Yd and Ya- 


( D / A )\ 114, 


0CC / X / \ _ 1 

P <M) W = E*i D, (r)(s <M) )„ t! A \r). (7) 

kl V 7 kl 

The Hamiltonian coupling is not Hda, but it is generally reported as the coupling between the 
Lowdin orthogonalized Yd an d Ya ■ F° r only two states this takes the form, 

Vda = [hda -Sda H dD ^ HaA ^ ■ (8) 

Turning to the superexchange picture, the effective coupling, is a summation of the contribution 
given by the interaction between D and A and the interaction of D and A with all bridges states, 
namely: 


Vda(E) =Vda+VdbE*b{E)\ba, 


( 9 ) 



where the superscript T stands for transpose, G g(E) is the Green’s operator, defined as 

G B (E) = -(Y B -El B )-\ (10) 

As shown in equation Eq. (9), Vda is the coupling for the donor-acceptor transfer, which in the 
absence of bridge states (CT through vacuum), would be the only contribution to Vda{E). On the 
other hand, if bridge states are present, the contribution to Vda{E ) is given by the second addend 
in equation Eq. (9). Generally, E appearing above is the energy at which the tunneling event 
occurs (i.e. at the crossing seam of the Marcus parabolas). In our works, 67 E was chosen to be in 
between Ea and £g, and specifically to be En Ea . This choice is invoked by several works in the 
literature 30,117-119 where it is well known that there is a mild dependence of the coupling with the 
tunneling energy. 118 However, this equation holds when there is no resonance between D, A and 
the bridges states. 23,30 ’ 120 ~ 123 If near-degeneracies appear then the transport regime transitions to 
resonant tunneling or hopping. 

2.1.2 Distance dependence of the electronic coupling 

In this section, we discuss calculations of the coupling matrix element (Vda) of hole transfer from a 
donor to an acceptor molecule through the vacuum. This means that the initial state of hole transfer 
is the donor molecule (D), and the final state the acceptor molecule (A), and no intermediate 
bridge states are considered. Any reliable method for computing couplings should be able to 
reproduce high level calculations of CT coupling in small molecular dimers. For this purpose, we 
initially chose 23 biologically relevant ^-stacks, 67 in order to analyze the distance dependence of 
the coupling, separations of 3-20 A were considered, as result a total of 276 coupling calculations 
were ran. Overall, our couplings show a good agreement with previous computations (e.g. we 
reproduced the decay factors, /3, ^-stacked dimers separated by vacuum). 

When a test set for hole transfer couplings featuring high accuracy couplings became avail¬ 
able 37 we could compare systematically the FDE-ET couplings with the benchmark values. 103 
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Benchmark calculations were ran on a set of 15 of ^-stacked dimers. This study was rigorous, and 
tested the effect of the basis set size, nonadditive kinetic energy functionals (NAKE) and exchange- 
correlation functionals (XC) on the value of the computed couplings. The most important finding 
that resulted from the benchmark work resided in the fact that GGA functionals coupled with a 
medium sized basis set and the PW91k NAKE functional allow the FDE-ET method to yield reli¬ 
able electronic couplings as tested against high-level correlated wavefunction (MRCI+Q, NEVPT2 
and SCS-CC2) methods applied to the array of dimers. The PBE and PW91 functionals are found 
to be a good choice in each case considered with a MAX error lower than 50 meV and an overall 
MRUE of a little over 7% in both cases. 103 Statistically, we found that hole transfer couplings are 
relatively insensitive to the choice of NAKE functionals, while our analysis of the basis set depen¬ 
dence shows that QZ4P basis set is the most problematic, as it often biases the FDE convergence to 
nonphysical states at short intersubsystem separations - a problem already well documented in the 
FDE literature. 124,125 Finally, Table Table 1 compares the performance of FDE-ET for different 
levels of theory. The results for GGAs are in good agreement with the benchmark values, and in 
some cases they showed to be superior to hybrid and meta-GGA functionals, particularly PBE and 
PW91. B3LYP also stands out as another valuable choice. 

Generally, all functionals perform well in the FDE-ET coupling calculations making FDE-ET 
a method that is relatively insensitive to the XC and NAKE functional choice. 

Table 1: Mean statistical values for the best XC-functional choices. PW91k is the NAKE along 
this work. Reproduced with permission from reference 103 


Set 

MUE(MeV) 

MRUE(%) 

MAX(meV) 

PB E/PW 91 k/TZP 

15.3 

7.1 

49.6 

PW 91 /PW 91 k/TZP 

15.2 

7.1 

49.1 

B3LYP/PW91k/TZP 

18.1 

7.9 

58.5 

M06-2X/PW91 k/TZP 

18.0 

8.2 

54.9 
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2.1.3 Hole transfer in DNA oligomers 


In this section, we discuss an interesting application of FDE-ET to charge transfer in biosys- 
tems.The electronic coupling for hole transfer in a completely dry B-DNA structure of G(T),yG 
and G(A),yG was calculated. The structures considered lack water molecules, metal counterions 
and phosphate linker groups. The latter is because the applicability of FDE is restricted to non¬ 
covale ntly bound molecular fragments. Consequently, appropriate modifications to the B-DNA 
structure had to be made: we have removed the phosphate groups and capped the dangling bonds 
with hydrogen atoms at 1.09 A from the bonding atom. The resulting structure of the modified 
G(T):vG is depicted in Figure Figure 4. The largest system considered is the double strand with 
ribose groups and counts 308 atoms and 1322 electrons. In this study, the role of the environment 
on the CT in DNA is elucidated and analyzed on the basis of an all-electron computation. 

Regarding the energetics (site energies), an uneven stabilization of the bridge states compared 
to donor/acceptor states occurs in both type of oligomers, being this effect more pronounced in 
the G(T)a,G system than in the G(A),vG system. By inspection of the overall electrostatics of 
the interaction between G:C and T:A, 126 we notice that T has a strong permanent dipole pointing 
towards A, similarly to C:G. Instead, A has a much weaker dipole compared to C or T and thus 
upon contact of the GTG strand with the CAC strand the cytosines will stabilize much more the 
holes on Gs than the adenines can stabilize the holes on Ts, hence the tunneling wall increases 
from single strand to double strand. 

Regarding the couplings, when the magnitude of the through space and through bridge cou¬ 
plings are inspected, our calculations show that the effects of the ribose groups and the nucleobases 
in the counterstrand are opposite and different in magnitude depending on the oligomer size (see 
table Table 2). We conclude, however, that the effect of the counterstrand on the computed su¬ 
perexchange couplings completely overpowers any effect due to the presence of the ribose groups. 
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Table 2: Through-space and through-bridge electronic coplings and tunneling energy gaps for 
single and double strand G(T),vG B-DNA, including the effects of the backbone (sugars). A - is 
shown for values below O.OlmeV. Reproduced with permission from reference 67. 



Vda (meV) 

^bridge ( OlC V) 

Edb (eV) 

Eba (eV) 

Single Strand no Ribose 

GG 

78.13 




GTG 

0.76 

12.46 

0.71 

0.50 

G(T) 2 G 

0.01 

1.13 

0.79 

0.66 

G(T) 3 G 

- 

0.09 

0.79 

0.77 

Double Strand no Ribose 

GG 

92.6 




GTG 

0.65 

7.66 

0.93 

0.96 

G(T) 2 G 

0.01 

0.47 

1.11 

0.94 

G(T) 3 G 

- 

0.02 

0.99 

1.16 

Single Strand with Ribose 

GG 

71.38 




GTG 

0.18 

25.01 

0.43 

0.37 

G(T) 2 G 

0.02 

1.70 

0.58 

0.37 

G(T) 3 G 

- 

0.21 

0.41 

0.41 

Double Strand with Ribose 

GG 

91.07 




GTG 

0.02 

7.35 

0.62 

0.87 

G(T) 2 G 

0.02 

0.61 

0.93 

0.60 

G(T) 3 G 

- 

0.02 

0.50 

0.82 
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Figure 4: The dephosphorilated G(T),vG B-DNA oligomer employed in the hole transfer coupling 
calculations. As the figure depicts, the hole tunnels from the bottom guanine (in balls and sticks) to 
the top guanine. The tunneling wall is provided by a series of three thymines (red branch, labeled 
as “bridge”). The counterstrand, C(A)a?C, acts as a solvating environment (in yellow, labeled as 
“spectators”) and no hole is allowed to localize on it. 

2.2 Constrained Density Functional Theory Applied to Electron Transfer 
Simulations 

Alternatively to FDE-ET methodology, constrained DFT (CDFT hereafter), a DFT-based proce¬ 
dure that was initially proposed by Dederichs et al, 127 and later introduced by Van Voorhis and 
Wu 128 with the aim of applying it to charge transfer reactions. CDFT is an effective method for 
calculating diabatic states for electron transfer, it relies on the idea of seeking the ground state of a 
system subject to a constraint. This can be achieved by adding to the conventional KS Fangragian 
an additional term that accounts for the constraining external potential, this reads as: 40 


■Z'cdft [p] = Ehk [p] + / v ext (r)p (r)dr - p 


p(r)dr — N e 


+ V C 


0 ) c (r)p(r)dr — N c 


( 11 ) 


same as regular KS—DFT 


where V c is the Fagrange multiplier of the constraint, G) c (r) acts as the weight function that defines 
the constraint, typically a population analysis based on a real-space 25 partitioning (such as Becke 
pop. analysis). N c is the value of the constraint, and at self consistency it should satisfy the 
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following tautology: 


N c = I ® c (r)p(r)dr 


( 12 ) 


Having defined the constraint parameters, the energy of the system can be computed by solving 
the KS equation for the constrained system: 

+ 1 V^\ dr ' +Vxc{r)+Vc(0cir) ) m (r) = ( r ) ( 13 ) 

where we have emphasized the functional dependence of the orbitals and orbital energies to the 
CDFT Lagrange multiplier. Clearly, the integral in Eq. (Eq. (12)) is only satisfied when an appro¬ 
priate choice of V c is employed. The term v xc is the exchange-correlation potential and </), are 
the KS-orbitals. Note that p( r) = 2 |0i( r )| 2 f° r closed shell systems. Thus also the density is 

a functional of the CDFT Lagrange multiplier. To our knowledge, the CDFT algorithm can be 
found on NWChem, 129 Q-Chem, 73 CPMD, 130 PSI, 131 SIESTA, 132 and ADF. 133 Computing the 
electronic coupling on a diabatic basis can be carried out similarly to Eq.(Eq. (5)-Eq. (8)) or using 
a CDFT-specific prescription. 41 

An example is the long range charge transfer excited states of the zincbacteriochlorin-bacteriochlorin 
complex (ZnBC-BC), an important structure in photosynthetic process in bacteria, has been cal¬ 
culated on the basis of CDFT procedure. 41,128 In figure Figure 5, the excited states at different 
intersubsystem distances is depicted, where the last point of each curve represent the CT excitation 
energy of the linked complex. These energies are in good agreement with previous methodolo¬ 
gies, 134 and also demonstrates that by constraining CDFT ground state the excitations are more 
accurate than TDDFT energies (1.32-1.46 eV). 128 

Additionally, CDFT can generate states with partial charges, 128 this is of particular impor¬ 
tance, for example in metal-ligand CT processes, where the diabatic states can be generated by 
constraining the charge on the ligand and metal center. 

Recently, the CDFT implementation of CPMD was tested against high-level wave function 
methods in the computation of electronic couplings for hole and excess electron transfer. 37,135 
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Figure 5: charge-transfer state energies of ZnBC-BC as compared to its ground-state energy at 
5.84 Aseparation. Lower line, Zn + BC ~;upper line, Zn BC + . Taken from ref. 128 

CDFT was shown to be on average within 5.3% of the benchmark calculations if 50% HF exchange 
was introduced (the average deviation goes up to 38.7% if HF exchange is not used). 

2.3 Fragment Orbital DFT 

The fragment orbital DFT or FODFT is a computationally low-cost method to calculate electronic 
couplings. This is because the wavefunctions of each diabatic state are approximated by the fron- 
teer orbitals of the isolated donor/acceptor fragments. 136-138 The underlying approximations in 
FODFT are that (1) the interactions between donor and acceptor have not effect on the orbital 
shape, (2) the coupling component related to orbitals below the fronteer is neglected (e.g. frozen 
core). In FODFT, the wavefunctions can be described by a single determinant of A — 1 spin-orbitals 
(j), where N — Na +Nd he. the sum of the number of electrons of the neutral donor and acceptor. 
These determinants are built from the KS orbitals of the noninteracting isolated donor and acceptor 
fragments. 



11 0.12 0.13 0.14 0.15 0.16 0.17 0. 

1/R (Angstrom' 1 ) 


* K (^ ■ ■ ■ <D ~ '* l “ c ) 
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The Hamiltonian used to calculate the CT matrix elements is the KS-Hamiltonian. Namely, 

n d +n a -i 

i=l 

N d +N a - 1 

H§* = I C (15) 

1=1 

where are the one-particle KS-Hamiltonians for either the "a" diabat or the "b" diabat. One 
feature of these Hamiltonians is that they are state dependent, thus, they are made of the combina¬ 
tion of orbitals of donor and acceptor species at the given state. The transfer integral, or coupling 
between states, is calculated as: 


H a ,b=(Ya\H\\l/ b ) 

(i6) 

Where N above is the fronteer orbital for D or A. Recently, Kubas et al. 37 have shown the dif¬ 
ferences of two FODFT flavors in the calculation of the hole transfer coupling for the HAB11 
database. As we can see in Table Table 3, the implementation including Na + Nd orbitals in the 
KS Hamiltonian (indicated by 2 N in the table) as done in ADF 136 is more accurate than the im¬ 
plementation using one of the Hamiltonians in Eq.(Eq. (15)) (which is indicated by IN — 1 in the 
table). 

FODFT has been successfully applied to models of CT in molecular semiconductors 139,140 
and also for modeling CT in biosystems. In the following, we provide applications of FODFT to 
biological CT: such as the determination of the hole rates on DNA hairpins linked by stilbenedi- 
carboxamide (Figure Figure 6), and the elucidation of the electron transfer between two cofactors 
in the SO enzyme. 
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Table 3: Hqa (meV) calculated with various FODFT approaches for HAB11 dimers at intermolec- 
ular separation of 3.5 A. Reproduced with permission from table XI of Kubas et al. 37 


FODFT(2N-l) 

FODFT(2N) 

ADF(2N) 

FODFTB 

REF 

Ethylene 

367.7 

389.2 

388.4 

343.7 

519.2 

Acetylene 

316.9 

345.8 

345.3 

212.0 

460.7 

Cyclopropene 

418.8 

443.7 

439.4 

367.4 

536.6 

Cyclobutadiene 

323.3 

346.9 

345.6 

261.6 

462.7 

Cyclopentadiene 

343.3 

360.6 

358.7 

283.2 

465.8 

Furane 

315.6 

334.0 

333.7 

280.3 

440.3 

Pyrrole 

328.7 

347.8 

347.7 

286.2 

456.3 

Thiophene 

341.2 

357.8 

356.1 

264.8 

449.0 

Imidazole 

310.7 

328.9 

328.2 

277.5 

411.6 

Benzene 

342.4 

353.5 

354.1 

299.9 

435.2 

Phenol 

190.5 

211.3 

279.5 

231.4 

375.0 


2.3.1 Hole transfer rates on DNA hairpins 

The absolute rates were determined by using Marcus theory, in Eq.(Eq. (1)), where the electronic 
coupling was calculated according to FODFT and the superexchange regime (see section Sec¬ 
tion 2.1.1). Knowledge of the rates forward and backward (see table Table 4) enables one to 
determine the equilibrium constant K = k t /k- t and the free energy change —AG = — keTln^K ). 
Comparable results with experimental results 141 were obtained. 

2.3.2 The Curious Case of Sulfite Oxidase 

An interesting and elusive candidate for electron transfer studies is the Sulfite Oxidase protein. 143 
For this protein, theory predicts an electron transfer rate between the cofactors (a heme and a 
molybdenum complex) about two orders of magnitude lower than what is measured experimen¬ 
tally. 2 To address this issue, Beratan et. al., using the Pathways model, suggested that the donor 
and the acceptor are joined together by a flexible tether. 144 As the tether allows the two cofactors 
to come sufficiently close to each other, electron transfer occurs at the rate shown by experiment. 
A recent simulation of this mechanism was carried out so that the protein was taken out of equi- 
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Figure 6: DNA hairpins used to calculate hole transfer via superexchange mechanism between a 
single guanine (green) and a guanine doublet (red). Reproduced with permission from Senthilku- 
mar et al. (Figure 1 of the reference). 142 


librium and positioned in a new folded state featuring a much decreased cofactor distance (about 
10 A). However, recent pulsed electron paramagnetic resonance measurements 143 indicated that 
the distance between the cofactors is unchanged on average from the one available in the crystal 
structure (32 A). To approach this problem using FODFT, the crystal structure of the protein is 
obtained and only the chains of the protein between the two cofactors are considered. The chains 
of interest are broken into individual molecules and treated as separate bridges. 


The FODFT computations that we present here will be part of a more in-depth study in a fu¬ 
ture publication. 145 Two ingredients are available from the simulations, the site energies, and the 
couplings between the sites. The energies of the hole transfer pathway for the electron transfer 
between the iron and the molybdenum is presented in Figure Figure 7. With the aid of Koopman’s 
theorem, the HOMO energies computed with FODFT are taken here as a measure of the ionization 
potential of each site. The simulation was able to shed light on some very interesting aspects of the 
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Table 4: Equilibrium Constant, K, and Free Energy Change, —AG, for hole transport between the 
proximal G site and the Distal GG Doublet. Experimental data in parentheses were taken from ref 
141. Table reproduced wth permission from ref 142 


Sequence 

K = k r /k- t 

-AG(eV) 

2b 

3.5 

0.032 


(7.5) 

(0.052) 

3b 

39.8 

0.093 


(6.7) 

(0.049) 

4b 

20.0 

0.076 


(> 0.5) 

(-0.02) 

5b 

10.0 

0.058 


(> 3) 

(> 0.028) 


couplings and the energy landscapes. Looking at the energy landscape in Figure Figure 7, it is clear 
that a path with highest couplings can be directly drawn between the donor and the acceptor. The 
landscape also shows the possibility of hopping stations-molecules that exist between the donor 
and acceptor and are close to them in energy. These two aspects of the landscape alone hint the pos¬ 
sibility of a hole transfer occurring over the 32 A. However, proteins are very complex structures 
with many variables such as size, dynamics and environment. Therefore, providing a quantitative 
analysis of the kinetic constant would require incorporating unbiased molecular dynamics and a 
more comprehensive structure to further characterize the role of these hopping stations. 

2.4 Ultrafast computations of the electronic couplings: The AOM method 

Recently, an ultrafast method to calculate electronic couplings was developed by Blumberger and 
coworkers. 135 The analytic overlap method or AOM is a useful method if CT simulations need 
to be coupled with molecular dynamics, like in proteins 146 or in organic semiconductors. 147 This 
quest requires hundreds or maybe thousands of Hoa an d site energy calculations. AOM offers an 
interesting alternative for such simulations. As in FODFT, AOM assumes that CT is only mediated 
by two SOMO orbitals (fronteer orbitals, similarly to FODFT), which correspond to each fragment. 
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Figure 7: Energy landscape for the hole transfer in Sulfite Oxidase: the circles represents the 
position of the center of mass of each fragment with respect to the electron transfer vector co- 
ordinate(the distance between the center of mass of the HEME complex (orange sphere) and the 
MOCO complex(pink sphere)); the size of the sphere is 1/ (x 2 ) in which x is the difference between 
the HOMO energy of the fragment and the Fe (that is HOMO of Fragment - HOMO of Fe). Fe and 
Mo were give a size of 1 for scale. 


Then, small Slater type orbital basis for the valence states is generated. Thus the overlap integral 
is evaluated as follows: 


Sda = ('PcI'Fa) = MWa) « 

atoms atoms 

~S DA = £ £ c *pn,i c pn,j{PnAPn,j) 

ieD jeA 

AOM further assumes contributions only from p-orbitals, particularly in organic compounds with 
7T—conjugation, the p-orbital considered is that one perpendicular to the plane of 7T—conjugation. 

Correlation of the overlap Sda and the electronic coupling given by FODFT is shown in Figure 
Figure 8 for a set of dimers and their geometries. In this picture, a satisfactory linearity between 
these two parameters is witnessed, therefore reliable proportionality constants can be achieved. 
H D a = CSda is used in order to obtain the constant C that can be used to get couplings of similar 
compounds. This first version of this method shows transferability for homo-dimers. However, 
transferability when non equivalent donor and acceptor systems are considered needs to be ex¬ 
plored. Nevertheless, AOM’s speed makes it a very valuable option. 
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Figure 8: Correlation between electronic coupling matrix element from sFODFT and overlap be¬ 
tween SOMO orbitals of donor-acceptor fragments. Taken from ref. 135. 


2.5 Note on orthogonality 

When carrying out a large number of coupling calculations, one encounters all those low prob¬ 
ability situations in which a method fails. In the case of FDE-ET, we probed a large number 
of so-called “difficult cases”. Specifically, we faced two limitations of the FDE-ET method. If 
the diabats are orthogonal, or quasi orthogonal, numerical inaccuracies arise in the inversion of 
the transition overlap matrix in equation Eq. (7). This is not specific to FDE-ET, but is a problem 
shared by all those methods that assume the diabatic states to be nonorthogonal. 30,115,148-150 when 
they are orthogonal, some of the equations previously developed simply do not hold anymore. Yu 
et al. have applied equations similar to Eq.(Eq. (5)-Eq. (8)) and obtained a picture of the behavior 
of the electronic couplings in the photosynthetic reaction center, see Figure Figure 9. If we con- 
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Figure 9: Correlation between electronic coupling matrix element and overlap between diabatic 
states for electron transfer. Taken from ref. 148. 


centrate on the left-side panels, we notice that in some cases the electronic coupling is proportional 
to the coupling, but in other cases (see lower left panel) the coupling seems to behave somewhat 
erratically as a function of the diabatic overlap. To understand this, let us consider two distinct 
limiting cases: (1) orthogonality by symmetry considerations, and (2) spatial separation of the or¬ 
bitals. We found that the second case is the predominant, as the distance between donor-acceptor 
increases the diabatic overlap becomes increasingly small. In the asymptotic limit, 44 there is a 
linear relationship between the coupling, the diabatic energy difference, and the diabatic overlap. 
If the overlap is small due to case (2), the asymptotic formula is not expected to hold. This explains 
the apparently contraddictory results presented in Figure Figure 9. 

Regarding FDE-ET, both cases can be circumvented computationally by performing a singular 
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value decomposition of the overlap matrix and then invert only those values which are larger than a 
threshold (i.e. Penrose inversion). For DNA presented in section Section 2.1.3, the default inversion 
threshold of 10 3 was appropriate in most cases. 67 However, three systems stood out: AG, GA and 
TT nucleobase pairs. All the systems above showed erratic behavior of the computed couplings 
for some specific donor-acceptor distances, specifically 4.0 A for AG, 3.5 A and 8.0 A for GA 
and 9.0 A for TT. We found that at those distances, the near singularity of the overlap matrix due 
to symmetry considerations (case 1 above) was the source of the erratic behavior. To circumvent 
these numerical issues, a threshold of 10~ 2 was adopted in these cases. 

We thus conclude that although there is a formal relationship linking the diabatic overlap with 
the value of the coupling at large donor-acceptor distances, 44 generally assuming linearity in the 
coupling vs. overlap (as mentioned briefly in the previous section regarding the AOM method) can 
lead to large errors in the magnitude of the computed couplings. In the future, inspired by a recent 
work by Evangelista et al., 131 a more stable algorithm that invokes an orthogonalization first, and 
then the computation of the couplings will be developed in our group. 


2.6 A fully semiempirical method: Pathways 

Pathways 2,50 is a semiempirical model which it is designed to reproduce electron transfer rates be¬ 
tween cofactors in proteins. 50,151 In essence, Pathways includes the contributions to the electronic 
tunneling from a stepwise path covering all nonbonded interactions, as well as the bonded ones at 
the nearest neighbor level. Namely: 


\h da \ 2 



(17) 


where e, are the steps the charge need to make from donor to acceptor. For example, a hydrogen 
bond is one of such steps. The above product is maximized by searching all possible steps that 
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contribute to the tunneling. The coupling is further split into three kinds of interactions: 


I H DA \ 2 = A 2 | Y\ £ bondi} 



space 



-bond 


(k) 


( 18 ) 


Pathways can yield reliable predictions of the electronic couplings, where the CT process in pro¬ 
teins are mediated by the interactions of a single or multiple configurations that the protein can 
adopt. 50 Pathways has been successfully applied to a number of CT processes in protein envi¬ 
ronment. For instance, the electron transfer between the proteins cytochrome c2 (cytc2) and the 
photosynthetic reaction center (RC) 152 in order to determine the protein structural dependence of 
this CT reaction, also, to look at the impact of structural and conformational variations on the elec¬ 
tronic coupling between the proteins methylamine dehydrogenase and amicyanin from Paracoccus 
denitrificans. 153 


3 High-accuracy electronic couplings 

This section is devoted to describing those methods which are able to predict the electronic cou¬ 
plings accurately given a certain definition of the corresponding diabatic states. These methods 
start with a mathematical definition of diabatic states (usually a definition that involves localization 
of the electronic structure) such that the resulting states resemble the donor and acceptor states in 
the electron transfer reaction. Once this is achieved, an adiabatic-to-diabatic transformation matrix 
is generated which can be applied to the adiabatic Hamiltonian to result in the diabatic Hamilto¬ 
nian featuring the sought electronic couplings in the off-diagonal elements. Usually, an accurate, 
wave function based level of theory is used for computing the adiabatic states and Hamiltonian. 29 
Examples of such techniques are, the Generalized Mulliken-Hush method developed by Newton 
and Cave, 29,30,154 Boys and Edmiston-Ruedenberg localizations of Subotnik et al., ,154, 1 55 and 
fragment charge difference proposed by Voityuk and Rosch. 42,156 Their utility lies on the possibil¬ 
ity of a very accurate computation of the corresponding adiabatic states, as was done for the hole 
transfer on 7T-stack DNA nucleobases at a CASPT2 and CASSCF level of theory accomplished 
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by Voityuk et. al. 35 That computation has served as the benchmark reference for many recently 
developed methodologies. 49,99,137,157 

Taking this as a motivation, let us briefly introduce each of the above methods, followed by 
some examples in which these methods where employed. 

3.1 GMH method 

In the two-state model, the charge localized diabatic states are related with the adiabatic states 
through the formula: 

£ 2,1 = 1 -(ed + E a ± \J(E d — E a ) 2 + 4 \H D a\ 2 ^J (19) 

where E D / A are the energies of the donor and acceptor states respectively, E\ and E 2 are the energies 
of the adiabatic states, this means the energy of the ground state (E\) and the first excited state (E 2 ). 
E 1 and E 2 can be obtained with any quantum chemistry method, however when highly accurate 
wavefunctions methods are employed also the resulting couplings will be of high quality. We now 
distinguish two cases: a symmetric case, for instance homo-dimers, and the general asymmetric 
case. In the symmetric case, we have that Ed = E A , and thus the electronic coupling does not 
depend on the diabatization procedure. Namely, 


2\H da \=AE 12 


( 20 ) 


AE \2 is the difference on energy between the adiabatic ground state E\ and the first adiabatic 
excited state E 2 . 

For asymmetric cases the GMH method prescribes that the proper diabatic states are those that 
diagonalize the adiabatic dipole moment matrix. In the two-state problem this is calculated as 
follows: 


\H da \ 


_ \l~in\AE 12 _ 

\j (Mil ~ M22) 2 + 4/i 7 2 


( 21 ) 
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where jU,y = ( v P/|/i £ 7 ’| v P J ), with p E r being the dipole moment in the direction of the electron 
transfer. The power of GMH lies on the way one calculates the adibatic states. As we have seen 
through all sections is that the authors benchmark their own method by calculating proper diabatic 
couplings by using ab initio methods as multireference Cl 37 and CASPT2. 35 

As an example of this method, let us discuss the very first example given by Cave and Newton 
on their paper. 34 The system Zii2H20 + , the transfer of a hole is done over the Zn atoms. However, 
the water molecule, which is located at a fixed distance opposite to the Zn distance, causes an 
energy splitting of the Zn orbitals. Thus the electronic coupling is determined for the following 
diabatic states: 


Zn{ l S) + Zn + ( 2 S) Zn + ( 2 S) +Zn( l S) s -»• s' (22) 

Zn( 1 S)+Zn + ( 2 S) Zn+( 2 S) +Zn{ l S) p -► p' (23) 

Zn(}S) +Zn + ( 2 S) Zn + ( 2 S) +Zn{ l S) s p (24) 

Zn(}S) +Zn + ( 2 S) Zn + ( 2 S) + Zn{ l S) p s' (25) 


where the transition are between the orbitals of the diabatic states (prime correspond to the acceptor 
state) (s-s’), (p-p’), (s-p’) and (p-s’). CASSCF wavefunction method was used in all calculations. 
In Table Table 5 we collect the values for the different couplings and /3,v. There the analysis of 
the distance dependence of the coupling is carried out for several Rozn distances. Note that for an 
infinite Rozn distance, the couplings for (s-p’) and (p-s’) are equal. 

We refer the reader to other publications which have evaluated the GMH method in detail in 
regards to its suitability in modeling two-state as well as multi-state problems for both excitation 
energy transfer and electron transfer processes. 39,156,158 

3.2 Other adiabatic-to-diabatic transformation methods 

Inspired by GMH, the electronic coupling can generally be obtained by rotating the corresponding 
adiabtic states into a set of diabatic states. Thus, each diabatic state can be expressed as a linear 
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Table 5: (a) Electronic coupling elements vs distance ( rz n Zn ) for ZniH20 + with r Zn o = 2.05 A. 
(b) Electronic coupling elements vs distance (r Zn Zn) forZ«2#2<3 + with r Zn o = 3-05 A. Results are 
in milihartree, the /3 values were calculated on the range of 5-9 A. Taken from reference 29. 


fZnZni A) 

Hss' 

H pp’ 

H sp ' 

Hp S ’ 

4.0 

28.3 

(a) 

23.6 

50.4 

42.7 

5.0 

10.5 

13.0 

51.7 

22.3 

6.0 

3.73 

7.55 

41.1 

10.1 

7.0 

1.09 

4.23 

21.8 

4.08 

8.0 

0.340 

2.57 

13.8 

1.62 

9.0 

0.0958 

1.44 

7.36 

0.611 

P 

2.28 

1.11 

0.81 

1.71 

4.0 

29.7 

(b) 

34.4 

59.3 

41.7 

5.0 

7.95 

14.7 

38.5 

22.1 

6.0 

2.34 

7.83 

19.1 

9.44 

7.0 

0.698 

4.25 

9.56 

3.84 

8.0 

0.203 

2.27 

4.78 

1.51 

9.0 

0.0558 

1.16 

2.32 

0.574 

P 

2.49 

1.32 

1.32 

1.74 
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combination of rotated adiabatic states as: 154 


Nstat es 

|s>= £ I<26) 

i -1 

Under specific assumptions of the nature of the system-bath interaction (the following is valid 
for the condensed phase), the coupling can be estimated by constructing diabatic states based on 
Boys, Edminton-Ruedenberg (ER) or von Niessen-Edminton-Ruedenberg (VNER) localizations. 
In Boys diabatization, the bath exerts a linear electrostatic potential on the system, thus the rotation 
matrix can be found by minimizing the following localizing function: 24,39 

Nstat es 

fBoys{U) = fBoys{^) = 52 = | (“/'l/JEr |—j) ~ (“/| ^ET |“/) | • (27) 

*,7=1 

Boys localization was shown to be equivalent to GMH for CT reactions. 39 

ER diabatization, dictates that the bath exerts an electrostatic potential that responds linearly to 
the field generated by the molecular system (sum of donor and acceptor) system: 


Nstat es 

/er{U) = /e/?(S) = 52 
i=l ' 


dr\ / dr 2 ~ 


,-|p(r 2 )|S / )-(S J -|p(ri)|S j ) 
|ri — r 2 | 


(28) 


In VNER diabatization, the bath exerts an electrostatic potential that responds linearly to the 
field of the total system, but the interaction potential is a Dirac delta function: 

^stat es r 

fvNER{U) =f VN ER{Z) = ~ £ / ^ r «S/|p 2 ( r ) |S/) - (S/|p ( r ) |S;) 2 ). (29) 

i=l 7 


Just like GMH, once diabats are generated, the electronic coupling readily arises from the off- 
diagonal element of the electronic Hamiltonian and is equal to: 


H D a = (Z D \H el \Z A ) 


(30) 


For a series of bridge mediated excitation energy transfer experiments, where the donor is 
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(a) Detachment density for the first lo 
calized diabat. 



(b) Attachment density for the first lo¬ 
calized diabat. 



(c) Detachment density for the second (d) Attachment density for the second 
localized diabat. localized diabat. 

Figure 10: Attachment/detachment plots for the occupiedaLS virtual separated Boys localized di- 
abatic excited triplet states near the avoided crossing. The molecule here is donor-C/A-acccptor. 
Taken from ref. 154. 


benzaldehyde and the acceptor naphthalene, 159,160 the transfer rates and couplings were calculated 
using Boys diabatization. In Figure Figure 10 it is shown how well the orbitals are localized on 
either donor or acceptor edges in the various diabatic states. Although this computation does not 
concern a CT process, we want to stress the ability of this localization procedures in generating 
true diabatic states. 


3.3 Fragment charge difference 


Similarly to the GMH method, the fragment charge difference (FCD) method yields a donor to 
acceptor coupling: 42 


\H da \ = 


|A#!2|AEi2 


(31) 


(Aqu -Aq 2 2) 2 +4Aq 2 n 

where Aqi and Aq 2 are the donor-acceptor charges differences in the respective adiabatic states \jj\ 
and \j/ 2 . Aq \2 is the off diagonal term and is defined in a general form as Aq t j = qij{D) — qij(A), 
i.e. the difference of the populations of the transition charges. 
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Finally, when donor and acceptor are in resonance, i.e when E D = Ea or Aq = 0. Hda = 
\{Ei-E 2 ). 

FCD method and its simplfied from (SFCD) are compared against GMH in the calculation of 
Hda for two Watson-Crick pairs GC and AT. In Table Table 6, the couplings calculated from FCD 
are in good agreement with GMH, the SCFD is also quite reasonable. Because the energy gap 
between donor and acceptor are large, the charge is completely localized on purines (lowest IP). 
However, if an electric field F (water molecule for instance) is tuned on near the pairs, the energy 
gap is then reduced and the coupling strength is enhanced. Overall, FCD is another good alternative 
to compute accurate couplings. However, the computation of the charges and the transition charges 
is dependent on the specific population analysis chosen. To our knowledge, only the Mulliken 
population analysis was used so far (i.e. the transition charges are evaluated on the basis of the MO 
coefficients over the atomic orbital basis set). 

Table 6: Hole coupling matrix element Hda of nucleobases within a Watson-Crick pairs. Energies 
in eV, dipole moment matrix elements in Debye, charges in a.u. Taken from ref. 42. 



Basis set 

E 2 E\ 

Ed~ Ea 

GMH 

H da 

SFCD 

Hda 

FCD 

Hda 

GC 

6-31G* 

2.163 

2.159 

0.0569 

0.0663 

0.0547 


6-311++G** 

2.092 

2.086 

0.0679 

0.0760 

0.0705 

AT 

6-31G* 

1.505 

1.502 

0.0421 

0.0524 

0.0363 


6-311++G** 

1.462 

1.459 

0.0474 

0.0528 

0.0425 


4 Practical Aspects: A protocol for running FDE-ET calcula¬ 
tions 

In order to obtain the electronic coupling for a CT reaction using FDE-ET, three different single 
point (SP) calculations have to be performed. FDE-ET is available in ADF. 51 In Figures Figure 11, 
Figure 12 and Figure 13, the input files corresponding to the FDE-ET methodology are described. 
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First, a single point calculation for each isolated fragment present in the system is carried out. This 
gives the initial density and energy of each subsystem without any interaction between them. It is 
important to save the check point files (TAPE files in ADF), because they contain all fragment in¬ 
formation needed in the subsequent calculations. Following all SP jobs for each isolated fragment, 
an FDE calculation is performed by taking into the account the whole supramolecular structure. So 
that, we create a diabatic state for each of the present subsystems. This is done by placing a charge 
different from neutrality in each subsystem, see Figure Figure 12. In this manner, two different 
directories are made: one in which an FDE calculation is carried out with subsystem 1 positively 
charged and one where subsystem 2 has the positive charge. In both cases, the SCF converges on 
the basis of subsystem DFT, thus, a series of three freeze-and-thaw procedure are done for each 
subsystem in each diabat. 

Once both isolated and embedded densities are obtained from the FDE calculations, a electron- 
transfer job is run whose purpose is to compute Eq.(Eq. (5)-Eq. (8)). As in the FDE calculation, 
the information about the fragment is of paramount importance in this electrontransfer job. In 
figure Figure 13, the input file that calculates the diabatic energies and the electronic coupling be¬ 
tween them is showed, in pink we can see that the check point files (t21.emb.rho* in the figure) 
corresponding to the embedded fragments in each diabat is copied directly to the ET directory, 
where the electrontransfer calculation is done. These files are renamed as fragA*.t21 for those 
ones from the diabat A (positive charge on the donor fragment) and fragB*.t21 for those that come 
from diabat B (positive charge on the acceptor fragment). It is worth mentioning that special care 
has to be taken in the management of the file names. As it is illustrated in Figure Figure 13, the 
fragments are numerated as 1 and 2, that means that the charge, departs from fragAl while fragA2 
is neutral and arrive to fragB2 while fragBl becomes neutral. This is very important when the 
system is comprised of more than two fragments, and the charge is moving throughout all of them. 
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$ADFBIN/adf « eor Bash Options 

Title : Fragment no. 1; (isolated) Job Titl© 

eprint prints the variables that you want in the outputfile 

SFO NOEIG NOOVL NOORBPOP 
END 

NOPRINT BAS FUNCTIONS 

units units to work with 

length angstrom 
a ngl e cfegree 
END 

symmetry nosym no groups of symmetry 
xc XC functional 

GGA PBE 
END 

geometry Type of calculation: Single Point 

SP 

END 

scf SCF cycles 

iterations 90 
END 

basis basis sets to be use 

Type TZP 
0re M> ne 
END 

atoms Cartesian coordinates 


1 

C 

0.000000 

0.000000 

0.000000 

2 

C 

1.332000 

0.000000 

0.000000 

3 

H 

-0.574301 

0.000000 

-0.928785 

4 

H 

-0.574301 

0.000000 

0.928785 

5 

H 

1.906301 

0.000000 

0.928785 

6 

H 

1.906301 

0.000000 

-0.928785 


END 

END INPUT 
eor 

Bash options 

mv TAPE21 t21.iso.rhol 

Figure 11: Single point (SP) calculation input file for an isolated ethylene. 
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$ADFBIN/adf << eor Bash Options 

Title : Fragment no. 1; (polarized) Job Title 

EPRINT 

SFO NOEIG NOOVL NOORBPOP 
END 

NOPRINT BAS FUNCTIONS 

units units to work with 

length angstrom 
angle degree 
END 

symmetry nosym no groups of symmetry 
xc XC functional 

GGA PBE 
END 

geometry Type of calculation: Single Point 

SP 

END 

scf SCF cycles 

iterations 90 
END 

charge 11 Charge and #of unpair electrons 

UNRESTRICTED 

fragments Subsystems involve in the calculation 

rhol t21.iso.rhol 

rho2 t21.iso.rho2 type=fde Frozen Subsystem 
END 


ATOMS 

Cartesian coordinates 

1 C 

0.000000 

0.000000 

0.000000 

f=rhol 

20 

1.332000 

0.000000 

0.000000 

f=rhol 

3 H 

-0.574301 

0.000000 

-0.928785 

f=rhol 

4 H 

-0.574301 

0.000000 

0.928785 

f=rhol 

5 H 

1.906301 

0.000000 

0.928785 

f=rhol 

6 H 

1.906301 

0.000000 

-0.928785 

f=rhol 

7 C 

0.000000 

3.50000 

0.0000000 

f=rho2 

8 C 

1.332000 

3.50000 

0.0000000 

f=rho2 

9 H 

-0.574301 

3.500000 

0.928785 

f=rho2 

10H 

-0.574301 

3.50000 

-0.928785 

f=rho2 

11 H 

1.906301 

3.50000 

-0.928785 

f=rho2 

12 H 

1.906301 

3.50000 

0.928785 

f=rho2 


END 


ALLOW PARTI ALSU PERFR AGS 

FDE FDE options 

PW91k 
XCNADD PBE 
END 


END INPUT 

Bash options 

mv TAPE21 t21.emb.rhol 

Figure 12: Single point FDE calculation input file for ethylene dimer. Both fragments rhol and 
rho2 come from two SP calculations for each isolated fragment. 
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#!/bin/bash 


fetch="#fetch frag??.t21" 
fetoh2="scp ../FDE/diabatA/t21 .iso.rho? 
cp ../FDE/diabatA/t21 .emb.rhol fragAl.t21 
op ../FDE/diabatA/t21.emb.rho2 fragA2.t21 
cp ../FDE/diabatB/t21 .emb.rhol fragBl.t21 
cp ../FDE/diabatB/t21.emb.rho2 fragB2.t21 

$ fetch 
$fetch2 

$ADFBIN/adf << eor 
Title electrontransfer run Job Title 

EPRINT 

SFO NOEIG NOOVL NOORBPOP 
END 

NOPRINT BAS FUNCTIONS 

UNITS units to work with 

length angstrom 
angle degree 
END 

SYMMETRY NOSYM no groups of symmetry 

XC 

GOA PBE XC functional 

END 

G s E p OMETRY Type of calculation: Single Point 

END 

SC T n SCF cycles 

iterations 0 ' 

END 

integration 6 . 06.0 Numerical integration 
unrestt? i cted Charge and #of unpair electrons 
FRAG ™ E -, NTS ^ n Subsystems involve in the calculation 

rhol t21.iso.rhol 1 

rho2 t21.iso.rho2 
END 

ATOMS Cartesian coordinates 

1 C 0.000000 0.000000 0.000000 f=rhol 

2 C 1.332000 0.000000 0.000000 f=rhol 

3 H -0.574301 0.000000-0.928785 f=rhol 

4 H -0.574301 0.000000 0.928785 f=rhol 

5 H 1.906301 0.000000 0.928785 f=rhol 

6 H 1.9063010.000000-0.928785 f=rhol 

7 C 0.000000 3.500000 0.000000 f=rho2 

8 0 1.332000 3.500000 0.000000 f=rho2 

9 H -0.574301 3.500000 0.928785 f=rho2 

10 H -0.574301 3.500000-0.928785 f=rho2 

11 H 1.906301 3.500000-0.928785 f=rho2 

12 H 1.906301 3.500000 0.928785 f=rho2 

END 

ELECTRONTRANSFER ET options 

numfrag 2 
InvThr 1.0e-3 
END 

zlmfit Density fitting 

END 

END INPUT 

eor 

Figure 13: Single point ET calculation input file for ethylene dimer. In pink there some bash-shell 
options in order to copy the information for each diabatic state calculated before with FDE. 


Bash options 

Densities and energies 
of both Diabats 
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5 Conclusions and future directions 


To conclude, we have presented our (fairly subjective) view of what tools are available nowadays to 
compute electronic couplings for charge transfer processes. We have surveyed in detail the FDE- 
ET method simply because we are among the developers of this method. Other methods based 
on DFT, and those that are best suited for being coupled with wavefunction based methods have 
also been discussed. The discussion also touches on the strengths and limitations of the various 
methods. 

When discussing the practical aspect of a coupling calculation, one must expose completely 
the methodology. We have done so for the FDE-ET method, and provided the reader with a step- 
by-step protocol on how to run such computations. This is important also for outsiders (such as 
experimentalists) as they can appreciate the kind of effort the theoretitians have to put in computing 
quantities relevant for the interpretation of the experiments. 

We apologize in advance to those authors who have developed all those methods that we have 
omitted from this presentation. Admittedly, we provide here a subjective view of the field. 
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